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Abstract 


Phase retrieval from the measurements of the Fourier modulus is an important and 
difficult problem. Among all the approaches developed to solve the problem so far, 
the iterative transform algorithms are currently the most efficient. However these 
algorithms suffer from major drawbacks that limit their practical applications. In 
this thesis the direct method i.e., the iterative transform method as well as wavelet 
based method of phase retrieval is being presented with the aim to discuss the 
means to improve the performance of the algorithms leading to better phase retrieval 
rjuality. 


It is bi'ing t'stablished that wavelet based method of phase retrieval increases 
the computational ('fhciency by over 60% in many cases while at the same time 
obtaining reduction in residual reconstruction error (upto 20%) as compared to the 
direct method. Further the quality criterion of wavelets chosen for subspace signal 
decomposition has an influence over the phase retrieval error and quality to a varjdng 
degree is being established. 
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Chapter 1 


Introduction 


Image reconstruction given partial information about its spectral contents remains 
an important and difficult problem. Although there is a one to one correspondence 
between a signal and its Fourier transform there are many applications [2] in which 
some of the Fourier domain information is either degraded or missing. The examples 
include Imaging through a turbulent atmosphere, X-ray Crystallography, Seismolog- 
ical and Ocean Acoustics applications and Astronomical Imaging/Applications. In 
order to restore the signal, therefore, it is desirable to recover the missing spectral 
information. An important area of investigation in such applications is to recover 
the complete loss of either Phase or Magnitude information. Also in the absence 
of any underlying signal model or constraints and also any apriori knowledge of 
signal, the loss of either Phase or Magnitude information of a complex function is 
considered irreversible. 

However although phase or magnitude alone is not sufficient, in general, to 
uniquely specify a sequence, yet under a variety of constraints an approximate signal 
may be synthesised in such a way that the information in original signal is fairly 
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represented in the synthesised signal even though the synthesised signal may not be 
an accurate representation on a point to point basis. The ability to reconstruct a 
signal from only phase or magnitude information alone would be useful in all such 
applications. To illustrate, for example in the case of problems referred to as Blind 
Deconvolution , a desired signal is to be recovered from an observation which is the 
convolution of the desired signal with some unknown signal. Since little is usually 
known about either the desired signal or the distorting signal, deconvolution of the 
two signals is generally a difficult problem. However, in the special case in which 
the distorting signal is known to have a Zero phase characteristics, the phase of the 
desired signal is undistorted. Such a situation is at least approximately known to 
occur in long term exposure to atmospheric turbulence. In this case the phase of the 
observed signal is approximately the same as the phase of the desired signal and, 
therefore, it is of interest to consider signal reconstruction from phase information 
alone. 

In this thesis the reconstruction of a 2- Dimensional (2-D) sequence from 
its Fourier Transform magnitude has been investigated in the framework of Unires- 
olution (Single grid) and Multiresolution Wavelet based approaches. It is supposed 
that the Fourier Modulus 1 xo(m.^)| of a roo,l and non negative object xo{m,n) is 
given. The reconstruction problem is hence to find a real and nonnegative unique 
solution x(m,n) so that \X{{j,,v)\ = 1xo(m, v)\- Since 

DFT[x{m,n)] = = \x{l^, i^)\ 

with I/)! = lxo(Ai, J^)l) tfio reconstruction of the object is equivalent 

to reconstruction of the Fourier phase (j>(iJ., u) from the magnitude (and hence the 
name Phase Retrieval). 

Many algorithms have been developed in the past to solve this problem[2,4] 
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and currently the most preferred algorithms are the Iterative Transform Algorithm 
and its variants. However these algorithms suffer from three major disadvantages. 

(a) Stagnation due to attractions to some local minima which causes the 
estimate to be trapped and subsequent non convergence, (b) Intensive computations 
due to the repeated requirement of 2-D DFTs at each cycle /iterations, (c) Lack 
of good criterion for choosing the initial guess to give the algorithm a better and 
unbiased chance of convergence. 

The platform for the current investigations are hence the Iterative Trans- 
form --Vlgorithm (hereinafter referred to as Fienup algorithm), its reconstruction 
properties, error convergence and quality of reconstruction on different classes of 
images. The variant of Fienup Algorithm i.e the Hybrid Input-Output has also 
been clubbed into the basic Fienup Algorithm for escaping possible local stagnation 
basins. The multiresolution analysis of the same images have been carried out and 
the parameters examined by using different wavelet bases. 


1.1 Phase Retrieval From Partial Information 

In electron microscopy, wave front sensing, astronomy, crystallography, and in other 
fields one often wishes to recover phase, although only intensity measurements are 
made. One is usually interested in determining an object, /(x), which is related to 
its Fourier transform, F(u), by 

F(u) = \F{u)\exp[ir!;{u)] = F[f{x)] 

— J-M f{x)exp{—i2Tru • x)dx 

where x is an M-dimensional spatial coordinate, and « is an M- dimensional 
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sj^atial frequeiu-y coordinate. For the majority of imaging problems however M=2. 
In practice for sampled data and assuming square arrays, u = {ui,U 2 ) and x = 
(.)•] . when' , H 2 , .ri, and X 2 = 0,1,2, N — l.Then one uses the DFT 

N-l 

F(u) ~ ^ f {x)exp{—i2'Ku ■ xfN) 

a :=0 

and its inverse 

N-l 

f{x) = N~‘^ ^ f {u)exp{i2TTu ■ x/N) 

w=0 

which are cumi>uted using the FFT method. 


I'or tiu' prohh'iu of recovering phase from two intensity measurements, as 
in ('h'ctron micro.scopy and in wave, front sensing, 

f{x) = \f{x)\exp[i'n{x)] 

is comph'x valued, and one wishes to recover 'ip(u) or equivalently recover r]{x) from 
UH'asurenient of both \F{u)\ and |/(a:)|. For the problem of recovering phase from 
a singU' inlr'usity measurement, as in image recovery from interferometry data in 
iistronomy and from structure factors in crystallography, one wishes to recover 
ar equivalently recov('r f{x) giv('n a measurement of lF'(w)| and the constraint that 
f{.r) be real and nonnegativf', 

f{x) > 0 


A i)articularly successful approach to solving these problems is the use of 
Fi('nup and related algorithms. These algorithms involve iterative Fourier transfor- 
mation back and forth between the object and Fourier domains and application of 
th(' measunHl data and known constraints in each domain. 

In what follows a generalised Fienup algorithm, alternative Gradient search 
nnd.hods and iterative Fourier transform algorithms are described for the case of 
single intensity measurements. 
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1.1.1 Error-Reduction Algorithm 


The predecessor of Fienup Algorithm, the Gerchberg-Saxton algorithm was orig- 
inally invented in connection with the problem of reconstructing phase from two 
intensity measurements(given intensity constraints in each of the two domains). 
The algorithm consists of the following four steps : (1) Fourier transform an es- 
timate of the object; (2) Replace the modulus of the resulting computed Fourier 
transform with the measured modulus to form an estimate of the Fourier transform; 
(3) Inverse Fourier transform the estimate of the Fourier transform; and (4) Replace 
the modulus of the resulting computed image with the measured object modulus to 
form a new estimate of the object. In equations this is ,for the A:th iteration, 


Gk{u) 

= !Gfc(u)|ea:p{i<^ife(u)] 


Ckiu) 



9k{^) 

= exp[i0;fc(x)l 


9m{^) 

= \f{x)\exp\i$k+i{x)] 

= \f{x)\exp[i9'k{x)] 


where gk-,dk,G'k^ and <()k are estimates of f,r},F, and 'tp respectively. As 
enumerated above the Gerchberg Saxton algorithm can be used for any problem in 
which partial constraints are known in each of the two domains. One only transforms 
back and forth between the two domains, satisfying the constraints in one before 
returning to the other. 

For the most general problem the error reduction algorithm consists of the 
following four steps: (1) Fourier transform gk{x), an estimate of f{x); (2) Make the 
minimum changes in Gk(tt), the resulting computed Fourier transform , which allow 
it to satisfy the Fourier-domain constraints to form Gfc(u), an estimate of f(u); (3) 
Inverse Fourier transform Ckiu)] and (4) Make the minimum changes in g'k{x), the 
resulting computed image, which allow it to satisfy the object domain constraints 
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to form gk+i{x), a new estimate of the object. In particular for the problem of the 
single intensity measurement (as in astronomy) the first three steps are identical to 
the first three steps of the Gerchberg Saxton algorithm, and the fourth step is given 
by 


9k+iix) = < 


0 , 


r 6 7 


where 7 is the set of points at which g'kix) violates the object domain constraints, i.e., 
wherever g^ix) is negative or where it exceeds the known diameter of the object. 
The diameter estimate can be taken since it will be roughly half the diameter of 
the autocorrelation function, which is the inverse Fourier transform of the squared 
Fourier modulus. The iterations continue until the computed Fourier transform 
satisfies the Fourier-domain constraints or the computed image satisfies the object 
domain constraints ; then one has found a solution, a Fourier transform pair that 
satisfies all the constraints in both domains. 


1.1.2 Steepest Descent Method 

An alternative approach to solving the phase retrieval problem is to employ one of 
the gradient search methods. The Steepest Descent method is closely related to 
the Error Reduction algorithm for the problem of reconstructing phase from single 
intensity measurement. 

An example of how this gradient search method will be used for our problem 
follows. One can define say B = Ep, the squared error in Fourier domain given by 
the following equations as the error metric which one seeks to minimise by varying 
a set of parameters. 
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B = jv-25:ig(u)-g'(u)P 

. u 

Here the iV^ values of g{x), the estimate of f{x) are treated as iV^independent 
parameters. Starting at a given point, gk{x), in the N^- dimensional parameter 
space , one would reduce the error by computing the peirtial derivatives of B with 
respect to each of the points g{x) {N^ partial derivatives ) forming the gradient 
and then move from gk{x) in a direction opposite to that of the gradient to a new 
point gk{x). One would then form a new estimate, ^^4.1(0:), from gl.{x) by forcing 
the object domain constraints to be satisfied . This would be done iteratively with 
suitable step sizes until a minimum is found (problems of minima whether Global 
or Local needs to be taken care of still). The above means that one minimises the 
error ,5, as a function of the parameters, g{x), subject to the object domain 
constraints. 

Ordinarily the computation of the iV^ partial derivatives would be very 
lengthy task since each evaluation of B involves an iV x iV DFT. However, for the 
problem of single intensity measurement considered here the computations can be 
greatly reduced by existing methods [6]. In effect the computation of g{x)t becomes 
more or less identical to the first three steps of the Error Reduction algorithm. 

The steepest descent method has been known to take more computations 
jtnd time for convergence than the Error Reduction algorithm. * 


1.1.3 Other Gradient Search Methods 

As stated in previous section for the phase problem of a single intensity measurement 
the Steepest Descent method is equivalent to the Error Reduction algorithm. Also 
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in practice the Error Reduction algorithm converges very slowly for the problem of 
single intensity measurement. Slower convergence of the Steepest Descent method 
has been observed for other applications as well [6]. Therefore the need of other 
gradient search methods. The Steepest descent method moves from the point gki^i) 
in parameter space to the point 

gl{x) = gk{x) - hkdgBk (1.1) 

where hk, the step size, is a positive constant, and the gradient is given by 

dgB = 2[g{x) - g (x)] 

where dgB is the partial derivative of error B with respect to a value at a given 
point, g{x). For many applications one would search along the direction of —dgBk, 
evaluating B repeatedly until the value of hk that minimises B is found; then from 
that point one would recompute the gradient and go off in a new direction. 

A useful block diagram representation of the family of gradient search meth- 
ods is shown as fig. l.l. 

A graduuit search method superior to the steepest Descent method is the 
conjugate-gradient method. For this method Eqn.1.1 is replaced by 

9k{^) = gk{x) + hkDk{x) 

where the direction Dk{x) is given by 

Dk{x) = g’kix) - gk{x) + {Bk/Bk-i)Dk-i{x) 

where one would start the first iteration with Di{x) = g[{x)-gi{x). After employing 
the above eqns., one would employ the spatial domain support constraints to form 
the new estimate as indicated in fig. 
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Figure 1.1: Block diagram- Gradient Search method 


A further method for this problem apart from the numerous other variations 
of the gradient search methods is the Newton-Raphson or the Damped Least Squares 
method. However each iteration of this method involves inversion of an x 
matrix while in the same time a large number of iterations of the gradient search 
method can be performed and therefore the Newton Raphson method although more 
sophisticated yet is not preferred. 

1.1.4 Input-Output Algorithms 

A solution to the problem of the slow convergence of the error-reduction algorithm 
has been the Input-Output algorithm, which has proved to converge faster for both 
the problem of two intensity as well as one intensity measurement. The input-output 
algorithm differs from error-reduction algorithm only in the object domain operation. 
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INPUT 


cr 

tD 


FMOD 


OUTPUT 


a 

o 


SATISFY 

FOURIER 

CONSTRAINTS 


-1 

FMOD 


Figure 1.2: Block diagram: Input Output concept 

The first three operations -Fourier transforming g{x), satisfying the Fourier-domain 
constraints, and inverse Fourier transforming the results- are the same for both 
algorithms. If grouped together as indicated in Fig. 1.2 , these three operations can 
be thought of as a nonlinear system having an input g(x) and an output g (x). The 
notable property of this system is that its output is always an image having a Fourier 
transform that satisfies the Fourier domain constraints. Therefore if the output also 
.satisfies the object <lomain constraints than that is the solution to the problem. 

It is evident that small change of the input results in a change of the output 
in the same general direction as the change of the input. More precisely, for a small 
change of the input the value of the corresponding change of the output is a constant 
a times the change of the input. Due to additional non linearity the corresponding 
prediction of output cannot be done exactly yet by appropriate change in the input, 
the output can be pushed in the general direction desired. If a change A 5 (a:) is 
desired in the output, a logical choice of the change of the input to achieve that 
change of the output would be Pl\gix), where is a constant (ideally equal to 
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For the case of Input-Output algorithm hence the next input is taken as(one 
of the many other choices) 

gk+i{x)= g'kix) + PAgkix) 

_ ^ gk{x), x ^7 

, xe'Y 

The use of above Eqn. is referred to as the Input-Output algorithm. If 
/? = 1 the same becomes the Error-reduction algorithm. 


Still another method of choosing the next input is 

I' 

9k+i{x) = < 

^ 9k{x) - .39k{^)^ ^€7 

The above is referred to as the Hybrid Input-Output algorithm which is an im- 
provement over all the above to avoid the stagnation problem. For this algorithm, 
if at a given value of x the output remains negative for more than an iteration, the 
corr(\sponding point in the input continues to grow larger and larger until eventually 
that output value must go to non-negative. 


1.1.5 Additional Methods 

Apart from the methods listed above the following other approaches are being used 
for Image restoration from incomplete data 

Restoration by Convex Projections 
Using Linear programming 
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Bayesian methods for estimation 


1.2 OVERVIEW 

The organisation of this thesis has been made as under 

Chapter 1 -Introduction, nature and scope of the problem. Existing algorithms for 
Phase retrieval from partial information. Alternative Gradient search methods 
and Iterative Fourier Transform methods are described. 

Chapter 2 -The Multiresolution Analyis 

Chapter 3 - Phase Retrieval 

Chapter 4 -Typical Reconstruction experiments, some practical considerations 
and results. 

Chapter 5 Conclusions and scope for future work. 
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Chapter 2 


The Multiresolution Analysis 

An important problem of image processing is to choose a representation that is 
well adapted for analysing the information content of the images as per the appli- 
cations. The sharp variations of signal amplitudes are generally among the most 
meaningful features. The information extraction and analysis is done more often by 
the decomposition of the image using various transforms. The earliest known such 
transforms were the Fourier representation. They use the orthogonal properties of 
the exponential family of basis functions and decompose the signal in terms of har- 
monic sinusoids. However they are characterised by their uniformity of scales in the 
time-frequency plane resulting in equal resolution analysis of the signal. When the 
signal includes structures that justify more close look into certain components of the 
signal, it is helpful to study the functions which have the microscopic ability and 
adaptability to say characterise the signals as per their frequency components on 
different scales. While short time Fourier transform STFT to a certain extent does 
the same signal analysis on a rectangular grid, the sampling of the time-frequency 
grid with wavelet transforms are very different. Lower frequencies, where the band- 
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width is narrow (the basis functions are stretched in time) are sampled with larger 
time stops while the higher frequencies are sampled more often resulting in a con- 
stant relative bandwidth as opposed to STFT. Hence the wavelet transform is a 
tool that cuts up data or functions into different frequency components, and then 
studies each component with a resolution matched to its scale. The wavelets had 
hence been constructed with a view ideally suited for multiresolution analysis and 
decomposition. 


2.1 Wavelets : A Review 

Wavelet transform of a signal implies the decomposition of a signal with a family of 
orthonormal bases m, n G Z, obtained through translation and dilation of 

a kernel function '>p{x) known as the mother wavelet i.e., 

- n) (2-1) 

Du(' to orthononuality, the wavelet coefficients of a signal / (.r) can be computed by 
f.lu' analysis formula 

C-rn,n = f (2-2) 

X 

To recover the signal from the wavelet coefficients, the synthesis formula is 

(2-3) 

m,n 


2.1.1 Wavelets : Construction 

To construct the mother wavelet ipix), the scaling function (j){x) with multiresolution 
analysis is used. The spaces of square integrable functions L^{R) is decomposed into 
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closed siibspaces l with the following properties 


• I 2 ’" C : The coarser subspaces are contained in the finer one. 

• n,n 6 j 2 l 2 "‘ = 0 : Separation condition. 

• U,ne 2 ’^ 2 "' = L'^{R)'- condition for completeness. 

• f{x) e V 2 -n /(2a:) e V 2 m+i : Scaling property. 

• There exists a scaling function (j){x) G V^o such that 

- n); m,n e Z (2.4) 

• Since <p{2x) G 1-21 C 1^20 and ^{2x) is a basis for the subspace 1 ^ 21 , the scaling 
function can be written as a linear combination of (f){2x) and it satisfies the 
two scale difference eqn 

4{x) = \/2j2h{k)(j){2x - k) (2.5) 

k 

'The wavehd, k<n-nel ‘tp{x) is given from the scaling function by the equation 

ij}{x) = ^/xJ 29 {k) 4 >{ 2 x - k) (2.6) 

k 

where 

g(k) = {-l)‘ft(l - k) (2.7) 

2.1.2 Wavelet Decomposition 


To perform wavelet decomposition, a recursive algorithm is implemented. The ex- 
plicit forms of </(x) and ip{x) are not required as they depend on h{k) and g{k) 
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respectively, given by Eqns., 2.5 and 2.6. A J level decomposition can be written 

<iS 


/o(^) = 

k 


J 

= + ( 2 . 8 ) 

k j-0 

where coefficients co,* are given and coefficients Cj+i,„ and dj+i „ at scale j + 1 are 
related to the coefficients Cj^k a-t scale j by the equations 


Cj^x,n = Y.Cj,kh{k-2n) ( 2 . 9 ) 

k 


and 


^j+i,n = Y1 ^j,k9{k - 2n) (2.10) 

k 

where O < J. thus equations 2.9 and 2.10 provide a recursive algorithm for wavelet 
decomposition through h{k) and g{k) and the final outputs include a set of J level 
coefficients !<;<■/ and the coefficient cj,n for the low resolution component 
<i>j,k(-^)- Similarly the recursive formula for the synthesis of the function is based on 
the wavelet coefficients dj^„; 1 < J < J and Cy,„ such that 


^j,k ~ ^ ^j+i,nk(k 2n) + ^dj^i^„9(k 2n) (2-11) 

n 


Fig. 2.3 illustrates the decomposition given in eqns 2.9 and 2.10. The pair 
of filters H and G with impulse responses h{n) = h{-n), g{n) = g[~n) are the 
QMF filters, corresponding to the half band low pass and high pass filters. Fig.2.4 
illustrates the reconstruction procedure implemented by upsampling the subsignals 
Cj+i and dj+i (interpolating) and filtering with h{n) and g{n) respectively, and sub- 
sequent addition of two signals. The wavelet transform decomposition is performed 
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recursively at the output of the low pass filter ‘H’. Thus the wavelet transform de- 
composes s signal into a set of frequency channels that have narrower bandwidths 
in the lower frequency region. 


2.1.3 Two Dimensional Wavelet Decomposition 

The wavelet decomposition can be likewise generalised to any dimension n > 0. 
Extending it to two dimensions for image processing, the signal is now a finite en- 
ergy function f{x,y) £ Let be the multiresolution approximation 

of The approximation of the signal f[x,y) at a resolution 2^ is equal to 

its projection on the vector space . Let be a separable multiresolution ap- 
proximation of Let 4>{x,y) — 4>{x)(j){y) be the associated separable two 

dimensional scaling function. Let be the one-dimensional wavelet associated 

with the scaling function cf>{x). Then the three wavelets are defined as 

0^(x,y) = </>(a:)v(y) 

0^(x,y) = ■tp{x)(f>{y) 

ip^{x,y) = 0(x)v(y) (2.12) 

Just as in the one dimensional case, the wavelet decomposition of a two dimensional 
signal involves convolving the signal with the separable, 2-D scaling function and 
the three wavelet functions given by equations 2.11. Just as the exact form of the 
scaling and wavelet functions in the one dimensional case was not required for the 
decomposition and the QMFs with impulse response h{n) and y(n) were utilised, 
the 2-D QMF filter coefficients are generated as follows 

hiL{k,l) = hik)h{l) 
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hHL{k,l) = g{k)h{l) 

hHH{kJ) = 9{k)g{l) 


(2.13) 


The 2-D wavelet decomposition and reconstruction are illustrated in the Figures 2.5 
and 2.6 respectively. 


2.2 Multiresolution Analysis: The basic idea 

Let A-zjfix) be the operator which approximates a signal at a resolution 2^. We 
suppose that our original signal f{x) is measurable and has finite energy ‘-fix) € 
L^(jR). Here we characterize A^j from the intuitive properties that one would expect 
from such an approximation operator: 

1. A23 is a linear operator If Azifix) is the approximation operator of some 
function f{x) at the resolution 2L then Azjfix) is not modified if we approximate it 
again at the resolution 2L This principle shows that Azi o A^j = .42;. The operator 
/l>;, is thus a projection operator on a particular vector space C L‘^{R). The 
vector space l'%j can be interpreted as the set of all possible approximations at the 
resolution 2-' of functions in L^{R). 

2. Among all the approximated functions at the resolutions 2^, Azi fix) is 
the function which is most similar to f (x) 

V^(x) € Il5(^) - ^ - /(^)ll (2-14) 

Hence, the operator Azi is an orthogonal projection on the vector space Vzj. 

3. Containment The approximation of a signal at a resolution 2^+^ contains 
all the necessary information to compute the same signal at a smaller resolution. 
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This is causality property. Since A 2 i is a projection operator on this principle is 
equivalent to 

Vj € Vi C ^2,+! (2.15) 

4. An approximation operation is similar at all resolutions The spaces ap- 
proximated should be thus derived from one another by scaling each approximated 
operator by the ratio of their resolution values 

Vjf G Z^f(^x) € Vi •<=>■ f{2x) G Vi+i (2.16) 

5. The approximation Azi of a signal f[x) can be characterized by 2^ samples 
per length unit When f{x) is translated by a length proportional to 2~^, Azifix) is 
translated by the same amount and is characterized by the same samples which have 
been translated. As a consequence of para 3, it is sufficient to express this principle 
for the resolution j = 0. The mathematical translations consists of the following 

• Discrete characterization. There exists an isomorphism I from Vj onto P{Z). 

• Translation of the approximation 

\fk € Z,Aifk{x) = Aif{x — k), where fk{x) = f{x - k) (2.17) 

• Translation of samples 

I{AMx))) = {ai)iez ^ (2.18) 

6. When computing an approximation of f{x) at resolution 2^ some infor- 
mation about f{x) is lost. However as the resolution increases to -t-oo the approx- 
imated signal should converge to the original signal. Conversely as the resolution 
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decrease to zero, the approximated signal contains less and less information and 
converges to zero. 


Since the approximated signal at resolution 2 ^ is equal to the orthogonal 
projection on a space V 2 > , this principle can be written as 


lim = 

J-4+0O 



is dense in L^{R) 


(2.19) 


and 


+00 

lim^V2,= f| ^2^ (2.20) 

j=-oo 

We call any set of vector spaces which satisfies the above properties, the multireso- 
lution approximation of L^(H). 


It has been seen that the approximation operator A 23 is an orthogonal pro- 
jection on the vector space . In order to numerically characterize this operator , 
we must find an orthonormal basis of Theorem 3.1 shows that such an orthog- 
onal basis can be defined by dilating and translating a unique function <;5(x). 

Theorem 3.1 Let be a multiresolution approximation of L'^{R). 

There exists a uniciue function q{x) € L‘^{R), called scaling function, such that if 
we set ^ 2 } (^) = 2 ^(j){ 2 ^x) for j £Z (the dilation of ^(x) by 2^), then 

( 2 , 21 ) 

is an orthonormal basis of V^ . For proof of above ref. [3]. The theorem shows that an 
orthonormal basis of any can be built by dilating a function <f>{x) with a coefficient 
2-' and translating the resulting function on a grid whose interval is proportional to 
2~K The functions <p 2 ^{x) are normalized with respect to the V-{R) norm. The 
coefficient \/2“-' appears in the basis set in order to normalize the functions in 
the L'^{R) norm. For a given multiresolution approximation {V 23 )jez, there exists 
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a unique scaling function (t>{x) which satisfies the above equation. However, for 
different multiresolution approximation, scaling functions are different. 

The orthogonal projection on can now be computed by decomposing 
the signal /(x) on the orthonormal basis given by Theorem 3.1. That is, 

-foo 

V/(x) € L^{R),A 2 jf{x) = 2~^ ^ - 2~^n) (2.22) 

n=— oo 

The algorithm for determining discrete approximation to the signal f{x) at 
resolution 2^ can be found as Mallat’s decomposition algorithm in [3]. Since (f){x) is 
a low pjiss filter, this discrete signal can be interpreted as low pass filtering of / (x) 
followed by uniform sampling at the rate 2T The scaling function (f>{x) forms a very 
special low pass filter since the family of functions ^v^ 2 “J 02 j(x — is an 

orthonormal family. 

Before we consider the wavelet representation, it will be pertinent to list 
one more useful theorem its under (for proof ref [3] ) 


Theorem 3.2 Let (p{x) be a scaling function, and let JT be a discrete filter 
with impulse response h{n) = Let H{uj) be the Fourier transform 

defined iis 


+00 

iy(a;)= X) 

n=-oo 

(2.23) 

H{uj) satisfies the following two properties; 


|JT(0)| = 1 and h{n) = 0(n"^) at oo 

(2.24) 

lJT(a;)| + |H(cx; + 7r)|2 = l 

(2.25) 

Conversely let H{uj) be a Fourier transform satisfying above equations and such that 

\H{co)\^0 for oj e [0, 7r/2] 

(2.26) 
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The function defined by 


0(a;) = n~ ,i/(2-Pa;) (2.27) 

is the Fourier transform of the scaling function. 

Our aim here is to build a multiresolution representation based on the dif- 
ference of information available at two successive resolutions 2^'^^ and 2^ . This 
difference of information is called the detail signal at the resolution 2K The ap>- 
proximation at the resolution 2-'+^ and 2^ of the signal are respectively equal to its 
orthogonal projection on V 2 J +1 and V 2 J. By applying the projection theorem, we 
can eiisily show that the detail signal at the resolution 2^ is given by the orthogonal 
projection of the original signal on the orthogonal complement of V 2 J iii ^ 2^+1 • Let 
O 2 , be this complement, i.e., O 2 J is orthogonal to K^i, 

O 2 J © V 2 } = V 2 J +1 . (2.28) 

To compute the orthogonal projection of a function f{x) on O 2 J, we need to find 
an orthonormal basis of O 2 J . Next theorem shows that such a basis can be built by 
scaling and translating a function 0(x). 

Theorem 3.3 Let be a multiresolution vector space sequence, 

<^{x) the scaling function, and H the corresponding conjugate filter. Let tp{x) be a 
function whose Fourier transform is given by 

(I) (I) ; G(u,) = (a, + tt) (2.29) 

Let 1^21 (x) = 2^ijj{2^x) denote the dilation of by 2^ then ^'\/2^'02i ~ 

is an orthonormal basis of and (y¥^’ip 2 i{x - orthonormal 

basis of i/ix) is called an orthogonal wavelet. 

We have hence seen that as (j>ik, K € Zis the orthonormal basis for V 2 j, any 
function belonging to L'^(R) can now be decomposed along all translates of d){2h). 
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The result is denoted by . Similarly the function can be decomposed along all 
translates of ip(2^x) and the result is denoted by D^. If the actual signal is having 
iV number of samples these results are having N/2^ number of samples. Thus the 
resulting signals belong to L'^(iY/2^). 

The decomposed signal along all the translates of (^(2-^x) is called the aver- 
aging signal and the decomposed signal along all the translates of ip{2^x) is called 
the detail signal. This is due to the nature of finite sequences h{n) and g{n) used 
to generate the <j){x) and ip{x) respectively. 


2.3 Quality Criterion For Wavelets Used In Image 
Processing 


Before choosing a wavelet for a particular application some essential quality criterion 
needs to be mentioned which characterise each type of wavelets [9] 

Wavelet regularity The regularity of the mother wavelet is important and appears 
to be closely related to the regularity of the signals to be processed. Since 
images are generally smooth to the eye, with the exception of occasional edges, 
it is appropriate to use regular wavelets. Indeed, there is a trade-off between 
wavelet regularity and the visual effect on the processed images. Because 
scaling function and the wavelet function ^(i) are related, they share the 
same regularity. Therefore in this discussion we restrict ourselves to the study 
of <p(t), which depends only on low pass filter H{z). Zeros at the Nyquist 
frequency play an important role in determining the regularity. One zero in 
H{z) at z = —1 is necessary to obtain continuity. To achieve regularity order 
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ol N, H{z) must have atleast N + 1 zeros at Z = —1. However if other zeros 
are also present, then they contribute towards killing the effect of regularity 
caused by zeros at ;; = — 1. This effect can be as much as 90%. This has been 
studied earlier and an algorithm to estimate the optimal regularity estimate 
is existing. The same is reproduced below: 

1. Remove all the zeros of H{z) aX z = —1. These zeros will contribute to a 
maximum Holder regularity of K, the actual number of zeros at 2 : = —1. 

2. The remaining factor F{z) = {1 + z~^)~^ H{z) will decrease this regular- 
ity. The decrease will be by the amount 1 — ct. 

3. (V, mentioned above is calculated by first iterating F{z) as under: 

F\z) = F{z)F{z^) . . . F{z'^^-^) 

Associated with this F‘^{z), we have time sequence /*. Then, oq is calcu- 
lated as under: 

CXi = 1/H0g2 maXn l/n+ 2 Jfc| 
k 

■1 Rx'gularity r is then given by 

r = K -IF ai (2.30) 

In Practice, i is iterated 20 times to get the regularity estimate. In general 
higher regularity of wavelets give better and tighter constraints. 

Vanishing moments It is the oscillatory character of the wavelet and is given in 
numbers N. For a wavelet ip{x) it is defined as 

= 0 yn = 0,l,...,N -1 (2.31) 

Since ^ is a wavelet, N is greater than or equal to 1 and, the property of the 
wavelets, where = 0 is verified. In fact, this property of number of 
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of vanishing moments is related to the regularity of the wavelet. That is, any 
r-regular multiresolution analysis of the wavelet will generate with r + 1 
vanishing moments. Also, all derivatives of upto the order N — 1 are 
zero at the point cj = 0. A wavelet with N vanishing moments enables the 
cancellation of all wavelet coefficients of a polynomial signal whose degree is 
less than *V. Thus, if / is a polynomial signal of degree less than iV, on support 
of i^rn,n then 

Cm,n(/) = (/, 0m, a) = 0 (2.32) 


Scaling function Spatial characterisation of the scaling function in terms of mo- 
ments which determine how 0(x) evolves with respect to x, which allows us 
to determine the energy concentration of 0 and provides information on the 
spatial length or the localisation of 0 is another significant aspect which needs 
mention. We introduce the criterion rrik which allows us to characterise the 
scaling function 0 spatially. To define m^t, we introduce the function p(x) such 
that 


p(x) 


i0(x)|' 


/ |0(x)12dx 

with p(x) > 0 and / p{x)d(x) = 1. ruk is then defined as 


(2.33) 


rrik 


y (x — mi)'^p{x)d{x) with mx = j xp(x) d{x) (2.34) 


Thus rrii is equal to the mathematical expectation of p{x) and m 2 corresponds 
to the spatial variance of the scaling function 0. m 2 allows us to determine 
the energy concentration of 0 and provides information on the spatial length 
or localisation of 0. This criterion also apply to the wavelet w{t). 


Characterisation of the associated filters To avoid distortion in image pro- 
cessing, the filter H{w) associated with the scaling function 0 must be linear 
phase or ideally zero phase. Indeed, non-linear phase filters degrade edges 
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and are more difficult to implement than linear phase filters. The number 
of elements making up the impulse response of h{n) must be small in order 
to limit the number of convolutions operations to be performed in the anal- 
ysis/reconstruction algorithm. It corresponds to wavelets whose support is 
compact (making the wavelet well localised). 
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Figure 2.1; Illustration of a Wavelet and their corresponding Scaling Function. Also 
illustrated is the Fourier modulus of the scaling function. Wavelet is a band pass 
filter. 
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Figure 2.2: Illustration of an orthogons 


responding Scaling Functions 


and a biorthogonal Wavelet and their cor- 





















Figure 2.5: 2-D Wavelet Decomposition 



Figure 2.6: 2-D Wavelet Reconstruction 
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Chapter 3 


Phase Retrieval 


3.1 Introduction 

Th(^ ivssue which is addressed in this chapter concerns the actual reconstruction of 
a multidinieusioual sequence from its Fourier modulus. It has been brought out 
in earli(!r chapters that in the absence of any apriori additional knowledge of the 
desirc'd seciuencc^ therc^ is as yet no practical algorithm which will always recover the 
correct phase from only magnitude information. Nevertheless, the Iterative Trans- 
form Algorithm has been discussed as also the Wavelet decomposition approach in 
solving the phase retrieval problem. Three different classes of images have been 
utilised for the same to investigate robustness of the procedures and results. 

Multiresolution analysis concerns itself with decomposing the signal spaces 
into subspaces using suitable wavelet bases, the quality criterion of which is to be 
considered for suitability in a certain application. Compactness, regularity and 
number of Vanishing moments play certain role in determining the efficacy of the 


32 



wavelet chosen for a certain task. To this effect many different length wavelet coeffi- 
cients have been used for multiresolution analysis. Also investigation for wavelets of 
three different chisses have been carried out namely Haar, Daubauchies and B-Spline 
Biorthogonal to see their effect on reconstruction quality. An effective tool during 
image compre.ssion and transmission would be the redundancy in the reconstructed 
image which have been checked in terms of truncation of Fourier coefficients of the 
reconstructed image. Low frequency component of the error responsible for the slow 
convergence of the iterative transform algorithm was another aspect which needed 
to be studied and therefore to broaden the picture it was seen as to the effect dif- 
ferent wavelet quality criterion have on this aspect and their effect on convergence 
and quality of reconstruction. 


3.2 Fourier synthesis from partial information: 
Importance of phase 


In this section wt‘ look at the importance of phase given an incomplete characterisa- 
tion of the Fourier domain. In particular, we will examine the the effect of combining 
the phase information with an estimate of magnitude function and the effect of com- 
bining magnitude information with an estimate of the phase information. Also the 
results obtained from course quantisation of the phase function will be investigated. 


It has been observed in a number of different contexts and applications 
that many of the important features of a signal are represented in the phase of 
its Fourier transform [2]. Specially, it is well known that a phase-only synthesis 
of a signal, formed by combining the phase of the Fourier transform of the signal 
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256x256 BIRO ORIGINAL 



(b) 


Figure 3.1; Phase only synthesis (a) original image and (b) the image synthesised 
using complete phase information and constant magnitude function 
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with either a constant or an ensemble averaged magnitude, contains a number of 
similarities to the true signal. From the reconstruction experiments performed on 
the test images reproduced in the figures[3.1]-[3.6] these are established conclusively. 
Shown in figure 3.1(a) for example is an original image and in Fig. 3.1(b) is the 
image formed by combining the phase of the Fourier transform of the image with a 
constant magnitude. 

y(m, n) = i/)}] 


On the contrast comparing this with a correct magnitude combined with 
zero Fourier phase(figure 3.4) image establishes that the essence of the image is pre- 
served in the phase of the image only.An even better phase-only synthesis may be 
obtained, however, by combining the phase of the Fourier transform of the image 
with a magnitude function which is more representative of images. 


Although much of the useful information in a signal appears to be rep- 
resented in the phase of its Fourier transform, the phase-only synthesis assumes 
perfect phase information. Yet it has been observed that even with a coarsely quan- 
tised tvpn'sentatiou of the phase, a reasonable synthesis may be obtained and the 
reconstruc:te<l image contains lot of similarities to the original as compared to the 
unknown phase.The same is being established with the help of test images (Figures 
3.2, 3.3). The one bit phase function is given by 


Q[(I>{H, u)] = < 


0 

TT 


if \<j>{u,v)\ < 7r/2 
otherwise 


As a result, the phase-only synthesis is simply given by 


where 


y{m,n) = v)] 


Y{u,v) 


1 i/l<^(/i,z/)| <7r/2 

— 1 otherwise 
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256x256 BIRD ORIGINAL 



(b) 


Figure 3.2: Magnitude only synthesis (a) original image and (b) the image synthe- 
sised using complete magnitude information and 1-bit phase function 
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256x256 CAMERAMAN ORIGINAL 
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An improved synthesis may be obtained by an ensemble averaged magnitude func- 
tion(as opposed to constant) which is obvious. 

Having observed the importance of phase information in the representation 
of images, it will be worthwhile to investigate the relative importance of magnitude 
information in signal representation. Therefore a set of experiments similar to those 
listed above were also performed to examine the importance of magnitude infor- 
mation. Shown in fig 3.4, for example , is an original image, x(m, n), and shown 
alongside is a magnitude-only synthesis formed by combining the magnitude of the 
Fourier transform of x{in, n) with zero phase 

y{m,n) = T[\X{ii,u)\] 

Clearly, there is no useful information portrayed with this synthesis. Also 
experimented was the combination of the magnitude of the Fourier transform of 
x(m, n) with the phase of another image (figures 3.5 and 3.6). As expected, the 
image produced by this .syuth(\sis bears a resemblance to the image whose Fourier 
phase wms used. 

The conclusion to be drawn from the above set of investigations are that 
there is, apparently, a significant amount of information conveyed in the phase of 
the Fourier transform of signals, which makes the phase retrieval problem very 
significant. 
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256x256 BIRO 0RK31NAL 



(b) 


Figure 3.4: Magnitude only synthesis (a) original image and (b) the image synthe- 
sised using complete magnitude information and zero phase function 


39 





(b) 


Figure 3.5: Magnitude only synthesis (a) original image and (b) the image synthe- 
sised using complete magnitude information of cameraman image and phase function 
from bird image 
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Figure 3.6: Magnitude only synthesis shows magnitude information of bird image 
and phase function from cameraman image 

3.3 Phase retrieval problem : The uniqueness ques- 
tion 

In this section we consider the reconstruction problem of a 2 - D signal in the back- 
ground of the uniqueness of reconstruction.lt has been mentioned earlier that with- 
out any additional constraints or information apart from magnitude information, 
a signal is not uniquely specified by the magnitude of its Fourier transform [4].For 
example, a signal may always be shifted (delayed) or multiplied by a factor of -1 
without changing the magnitude of its Fourier transform. Thus, any signal y(m, n) 
which is related to x(m, n) by 

j/(m, n) = ±x(m + k,n + l) 
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for any integers k and I has a Fourier transform with the same magnitude as x(m, n). 
In addition, rotating a signal by 180° [or reflection about the point (0, 0)] will not 
change the magnitude of its Fourier transform, so x[m,n) and 

y{m,n) = x{-m,-n) 

have the same Fourier magnitude. In most applications, perhaps these ambiguities 
are of no concern since it is the relative amplitudes of the various samples within 
the signal that are of primary importance. Therefore although such reconstructions 
have been observed in many reconstructions in this thesis, yet two signals x(m, n) 
and y{m, n) will be said to be equivalent as per convention if they are related to 
one another by any combination of the above relations.Excluding this ambiguities, 
it is still not possible to uniquely define a signal in terms of its spectral magnitude 
because of the following: 

• It is always possible to convolve s signal with an arbitrary all-pass signal(say 
one which has a Fourier transform with unit modulus also) to obtain another 
signal with the same spectral magnitude. 


3.3.1 Uniqueness :Revisited 

Having said, that without any additional information or constraints a signal is not 
generally defined by the magnitude of its Fourier transform, with only the finite 
length constraint, we state that due to the irreducibility of almost all polynomials 
in two or more variables, an important conclusion can be drawn form the unique- 
ness theorem [4], [2]. In particular it has been shown in the above theorem that a 
2-D signal with finite support which has an irreducible x-transform is uniquely de- 
fined(to within a sign, time shift, and time reversal) by the magnitude of its Fourier 
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transform. Although this theorem assumes that the Fourier magnitude is known for 
frequencies, yet this can be generalised to show that uniqueness still holds if the 
magnitude of the DFT of x(m, n) is known. 


3.4 Direct Approach:The Fienup Algorithm 

The image reconstruction problem (also known as phase retrieval problem) can be 
stated as follows: 

Given the magnitude of the Fourier transform,lF’(u, u)|, of a bandlimited image 
f{x,y) € L^{R), and some apriori information on /(x, y), reconstruct the image 
/(x, y), or equivalently, retrieve the phase $(u, v) of F{u, v). 

To this effect the Fienup algorithm was chosen due to its wide acceptability 
and practicability. It was supposed that the image Fourier modulus was the only 
information given and any other apriori information was assumed as not available. 

The Fienup algorithm is based on the idea that by iterating between the 
spatial and frequency domains and successively enforcing the known constraints at 
each stage one can obtain a numerical solution that minimizes the distance between 
the measurement and its estimate. For the most general problem the Fienup algo- 
rithm consists of the following four steps(for the k^/^ iteration) 

1. Fourier transform gkix,y), an estimate of /(x,y), yielding Gk{u,v) 

2. Make the minimum changes in Gk{u,v) which allow it to satisfy the Fourier- 
domain constraints to form GKu.v), an estimate of F(w, v) 

3. Inverse Fourier transform .yielding 9k{x,y), the corresponding image 
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4. Make the minimum changes in pjt(x, y) that allow it to satisfy the object 
domain constraints to form y), a new estimate f the object. 

For phase retrieval for a non negative object from the Fourier modulus 
|F(u, u)|, these four steps are 


1 . 

2. 

3. 

4. 


Gk{u,v) = \Gk{u, v)\exp[i(f)k{u,v)] = T[gk{x,y)] 
G'k{u,v) = lF(«,u)jexp[i<?!»fc(u,t;)] 

= -F~^[G]t(n, u)] 


9k+i{x,y) = < 


9ki^,y)r 

0 , 


(a:,y) € 7 


Where 7 is the set of points at which g^{x, y) violates the object-domain constraints 
and where y*, and and 0 * are estimates of /,F,and the phase ip of F, respectively. 


The algorithm was started by using an array of random numbers uniformly 
distributed between 0 ~ 255 assuming an 8 -bit format for further image processing. 
In this thesis as has been mentioned earlier the step four of the algorithm has been 
modified as per Hybrid Input-Output algorithm for getting better stagnation profile 
during reconstruction.Hence the step 4 has been taken as under: 


9k+i{x,y) = < 

k 


9k{x, y) - pg'kix, y), (x, y) € 7 


Where is a constant feedback parameter which controls the level of change in the 
next iteration updated guess. Values of ^ can vary between 0.5 to l.For this thesis 
the value of as 0.7 worked well and was hence fixed for all reconstructions. 
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In practice a number of cycles of iterations were performed, where one cycle 
consisted of, 40 — 50 iterations of the Hybrid Input-Output algorithm followed by 
5 - 10 iterations of the basic Fienup algorithm. 

For the general images including astronomical images (Available 8-bit 128 x 
128 data of Saturn taken by Hubble space telescope was used as a test case), the 
object domain constraints are the objects nonnegativity and a (usually loose) sup- 
port constraint. Then the points in the set 7 are those outside the assumed support 
and thOvSe within the support for which g'k{x^ y) < 0. The diameter of the object was 
estimated since it is just half the diameter of the autocorrelation. However it is evi- 
dent that such estimation is only rough and the exact support of the object cannot 
be determined uniquely from the support of the autocorrelation and so the support 
constraint cannot be applied tightly. At best the support can be subsequently used 
for cleaning up the image if satisfactory reconstruction has been obtained in the first 
run. Computational expense has to be however kept in mind for such inadequate 
support constraints. 


3.4.1 Convergence for the algorithm 

\ 


A measure of convergence of the algorithm to a solution (a Fourier transform pair 
satisfying all the constraints in both domains) is the normalised root-mean-squared 
(NRMS) error metric in the Fourier domain: 


Fp = 




1/2 


( 3 . 1 ) 


or in the object domain. 


Eo = 




1/2 


( 3 . 2 ) 
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where 7 is defined as earlier. The summations are performed over all points in 
image or Fourier space. Since the convergence of the algorithm with little known 
constraints is usually slow, it is only natural that error be noted only after atleast 
one cycle of iterations comprising of say total 40 - 100 iterations of the Fienup 
algorithm. 


3.4.2 Additional aspects of the Fienup algorithm 

Most often the Fienup algorithm will give a reasonable image reconstruction and 
when the error goes down to within limits of visual satisfaction, one can think that 
the solution has been found. Alternatively as the error starts stagnating for atleast 
3-4 cycles than one can assume that the reconstruction is trapped in a minima. It 
is unlikely that error will ever go to zero because (apart from others) 

1. The Fourier modulus resolution may not be satisfactory to resolve the recon- 
struction. 

2. Tight or incorrect support constraint makes the reconstruction inconsistent 
with either the nonnegativity constraint or the Fourier modulus. 

3. Choice of initial input to the algorithm varies the error profile as has been 
observed [4] 

For the general problem one has a nonnegativity constraint in the object 
domain. Furthermore one can compute the bounds on the support of the object 
from the autocorrelation which roughly resembles that of the object if the Fourier 
modulus resolution is sufficient. The simplest way is to use a rectangle that is half the 
size, in each dimension, of the smallest rectangle that encloses the autocorrelation. 
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Actual support of the object has been considered to be not known in our case and 
thresholding of autocorrelation to .005 of the maximum value has been used to find 
a support mask. The support mask has been defined as a mask which is unity inside 
the support and zero outside. In effect the autocorrelation is thresholded to .005, 
demagnified by a factor of 2 in each dimension by discarding every other row and 
column and filled up with ones to construct a mask. 

There are many ways to pick an initial estimate as an input to the edgo- 
rithm [2] and it has been established that although same crude phase estimation 
methods may have worked with some reconstruction problems but the efficacy has 
not been established for all. On the other hand it has been found that random 
numbers perform as well. In case some other known apriori information about the 
nature of phase is there then that should be used since having an initial input close 
to the true solution not only reduces the computational complexity but helps to 
reduce the inevitable stagnation problems. To give the algorithm an unbiased start 
in the absence of any other information, random numbers were therefore chosen as 
the initial input in the object domain(alternatively in the Fourier domain, random 
phase could have been chosen). 

The algorithm can be made to converge faster and avoid a stagnation prob- 
lem if the support mask is chosen to be somewhat smaller than the correct support 
for the first few cycles. Since it is the incorrect support, the smaller support mask is 
inconsistent with Fourier modulus and stagnation eventually occurs. Nevertheless, 
the smaller support mask helps to force most of the energy of the output, g {x,y), 
into a confined region in fewer iterations. After this has happened the support mask 
should be enlarged to the correct support constraint as given by the known or esti- 
mated value. During the last few cycles or alternatively if stagnation occurs{error is 
being computed at the end of each cycle), it is beneficial to enlarge the support mask 
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to larger values to give breathing space to the reconstructed image. In addition it 
helps since recall that the reconstruction is always having the trivial ambiguities(of 
translation and reversal) and therefore the support mask often starts hurting the 
reconstruction by the way of truncating the partially reconstructed image inadver- 
tently. Keeping all these factors in mind, the algorithm was tested with an adaptive 
mask which increased every four cycles. Satisfactory reconstruction error, stagnation 
profile and visual quality were hence obtained. 

The discrete Fourier transforms are computed with the FFT algorithms. 
The sampling in the Fourier domain should be fine enough to ensure that the ob- 
ject domain array size is atleast twice the dimensions of the object itself, which is 
equivalent to achieving the Nyquist sampling rate for |F(u,u)p. 

A practical method of implementation of the algorithm is to compute the 
phase from the real and imaginary parts of Gk{u,v), then combine it with |F(u,u)| 
to form G'kiu, u),and finally compute the real and imaginary parts of G'k{u, v) (which 
are required by the FFT) from its modulus and phase. 

Depending on the many factors listed above the Fienup algorithm converges 
to some solution and in some cases requires some additional inputs.The test case of 
‘Claire’ image was one such experiment which did not give a reasonable reconstruc- 
tion. Apparently the Fourier resolution was not adequate to resolve the details and 
hence the outline of the image was only visible along with some poor visibility inside 
the body of the image. It seems increasing the Fourier resolution only marginally 
helps the reconstruction and there are classes of image which are difficult to recon- 
struct. Additionally among others, one reason for such poor reconstruction could 
be the absence of the detail coefficients in the reconstruction which are responsible 
for blurring of the edges. 
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Some other problems of reconstruction and stagnation has been adequately 
documented in literatures [2], examples of which are 

• Appearance of simultaneous twin images in reconstruction. 

• Superimposed stripes on reconstruction. 

However these problems were not encountered in the test cases taken in this thesis. 


3.5 Wavelet method :The Multiresolution approach 


A major factor in the non convergence of the Fienup algorithm and other related 
iterative transform algorithms is the choice of the initial input to the algorithms. 
In general it is known that the closer the initial input is to the true image, the 
better will be the cotivergence and subsequent less chances of stagnation. Another 
major disadvantage in the iterative methods in general is the intensive computa- 
tional rt'cjuirements of the pair of 2-D DFTs required at each iterations of a cycle 
of the algorithm. Since the iterative transform method operates only on a fine grid, 
the computational complexity is enormous and any saving on this part should be 
investigated. The wavelet decomposition approach can significantly improve the 
performance of these algorithms in most of the cases in which reconstruction is gen- 
erally obtainable from the Direct methods. In so far as the visual reconstruction 
quality is concerned, it was found that a marginal improvement in error profile does 
not apparently improve the image quality; which may also be due to the Umitation 
of the eyes to discern. Nevertheless it is conclusive that more the input estimate is 
closer to the true image, the better should be the chances of an improved construc- 
tion. As has been mentioned earlier the quality of reconstruction, the error profile 


49 



aiul tha stagnation profile, all depend on the initial input in some way or other, the 
nit<'nTlat ionship lu'tweiai them is not clearly understood. This results in an even 
iiit’rc>as(‘<l error in a lew inst ances at some cycle of iterations vis-a-vis the class of the 
image. Ment ion of such instances has also been found in literatures [6]. Nevertheless 
tin' mull iresolut ion wavelet approach provides a rough(and quick) estimate of the 
solution at a low resolution that is later refined employing coarse to fine strategy. 
Folh>wing; advantag('s accru(^ out of this approach: 

1. Fnahles the algorithm to avoid stagnation by providing a better initial guess. 

2. improving t he convergence rate by decomposing the search space into orthogo- 
nal subspact's and therefore reduction of low frequency component of the error 
responsible for the slow convergence of the algorithm. 

3. RediK'tion in computational complexity since the number of variables to be 
{‘stimated at the coarser level is lesser than that at full resolution grid. 


3.5.1 The algorithm 

L(‘t fir) be a signal with finite support and let Am be the operator that approximates 
a signal !{x) at resolution m. Without loss of generality, we can assume that 
f{x) € Vy where is the linear spanned by the scaling functions (pm,n- Hence, we 
can oomput,e /"‘ € V),, the discrete approximation of f{x) at any resolution m [3] 

by 

= Am f = 

00 . s, 

= H‘^n-k)r-\k) (3.3) 
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whvrv is ih<> scaling function generated by the dilation equation 

OO 

'r'mti,»{-r) = i‘P^n+l,k{u),(prn,k{u))(fim,k{'^) (3.4) 

k=r.^OO 

and U is the mirror iilter (i.e. /r(n) = h{~n)) of the finite impulse response(FIR) 
filter h\ii] defined by 

^ r+oo 

h(n) = v2 / ifi{x)ip{2x — n)dx (3.5) 

7—00 

H> taking the inner product of / and (Pm+i,n(^) in eqn. 3.4, and rearranging the 
t(‘nns. (Uie can show that 


OO 

= E h{2n-k)r{k) (3.6) 

A;=— OO 

'Hk' alu»vt‘ in effect means that 


/'"+^(n) = h(2n) * /”^(2n) 


If a coarse approxiniation /"‘ € Vm of f{x,y) is obtained, it can be used as an 
initial estimate for the iterative algorithm to obtain a new finer approximation 
/"* ’ < r,„ i of This procedure is being repeated until m = 0, where the 

image is iH'constructed at the original resolution.In the experiments under taken, 
normally 2 3 lev<‘Ls of decomposition were under taken. It is to be understood 

ht‘re thiit the support constraints are just as valid for the decomposed subspaces 
as was for th(‘ original space of full resolution(m = 0). This in effect introduces a 
constraint of computing the subimage support from the subimage autocorrelation 
since t h«' suhspawi image dimensions are different from half the original image di- 
nnmsions.l'hen'fore there was a limitation in terms of how lower (coarser) one can 


go to get c;orr(‘ct solutions at the coarse levels.The alternative which was chosen 
was that the support at child subspace is taken approximately half the mother sub- 


space.Though this results in some truncation of the reconstruction at the subspace 
l(!V(d, yet to sonu^ extent an adaptive mask takes care of this. 
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:£^*TRAL library 

1. 1. T., 




i»*A mn^ 



In onlnr tn n'constnicl wo need to use the Fourier transform magnitude 
iU th*' rrsulutioti. Hinco a precise measurement of those are not available at 
all lit!' required n‘s>lutinns. i.o„ jF"*(w,u)j;0 <m < M,an estimate of these mea- 
siirruKuiis is obtained by taking the Fourier transform for both sides of eqn., 3.6 as 
foiiowh; 

( 3 - 7 ) 

neitee. iht' Fourier trati.sform magnitude at all resolutions m > 0 was ob- 

tained from liie Fourier transform magnitude of the image at the original resolution 

Hi f). 


I he ;!i.;.!H\i[uali.iu» /"* G were then obtained by performing several 
eyeles of the basic Fiemip algorithm at this resolution. The procedure was typically 
started by e-eoid'-i in;; an array of 16 x 16 random numbers(since resolution level 
Hi (K 1. 2 and 3 were considentd in most of the cases)as the initial guess as input to 
llie algtu'ithm. Tlie mitocorrelatiou was computed from the Fourier modulus at the 
original resolution and thn'sholdcd at 0.005 of the maximum value for compactness, 
else muse if jiny would have iulluenced the support region resulting in grossly wrong 
reetuisimetitius ntnl stagmition. An ad hoc value ol 0.005 was found to be giving a 
satisfactorv support hmgth with the constraints of the Fourier resolution. Once the 
estimat** <tf the image is obtained at the coarser resolution / , m > 0, the resulting 
image is expanded by a factor of two by interpolating rows and columns with zeros 
and preserving the dimensions of the next fine image subspace. The magnified 
version ’ beeom<‘s the initial input for the next higher resolution. Starting with 
this new initial estitnate and the Fourier magnitude at this resolution 
the same bmsic steps are performed at this resolution level m - 1. This coarse to 
finv. proec'SH is repeated until the image is reconstructed at the original resolution 

m = b. 
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In terms of computational efficiency, it is clear that if the number of pixels 
in tht' original n 'l ci iun image was N x N, then the approximate subspaces at 
each c'03irs.‘r resolution levels(^,„/ if original image was /) wiU have 2-^N x 2"”* AT 
pixels, biiue the iterative tr<msform algorithms are computationally intensive, it is 
a eonsideriible saving if one can work in the subspaces. 

The normaii.sed root mean square reconstruction error(NRMS)was hence 
calculated akin to the procedure as explained earlier(section 3.4.1) for all the wavelet 

cL'vsse.s. 


3.6 Ambiguity Correction 


As has been bought out in earlier sections(see 3.3), the reversal and translation are 
two ;t!!ibignifii s which are inherent in the phase reconstructed images when the only 
known inptit to start with is the Fourier modulus of the original image. This is due 
to the fa<’t that the Fourier modulus of an array is equal to the Fourier modulus of 
thtt time .shiftwl vt^rsum ;is well as origin flipped(reversed) version array. Since phase 
retrieval from the wjivelet decomposition method employs multiresolution technique, 
an attempt was made to correct the ambiguities before the phase retrieval actually 
starts at the full rc.solution on the finest grid. Hence at the coarsest level the re- 
construction was vi.simlly seen and the ambiguities corrected in the following manner: 


1. Reflection ambipity - Two signals a;(m,n) and y{m,n) will be said to 
be in reflection ambiguity if they are related by 2 /(m, n) = x{-m, -n). The Fourier 
modulus of both signals x{—m, — n) and y{m,n) are same but they are the phase 
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reversed version of each other which can be easily found out by the basic equation 

•fCX) 

J^[y{m,n)]= y(m, (3.8) 

m,n=--oo 

Therefore the reconstruction was visually seen at the coarsest level and and reversal 
ambiguity if present was removed by reversing the phase of the subsignal obtained at 
that point. The iterations were thereafter allowed to proceed further with the phase 
reversed signal. It is notable that the algorithm preserves the phase of the signal 
in whatever is reconstructed the first time without further reversing or translating it. 


2 Translation ambiguity - This is manifested between two signals in the 
form y(m, n) = x{m + k,n + 1). Again we note that the two signals have the same 
Fourier modulus. The phase of the two differ in a sense that translation in spatial 
domain is reflected as a modulation term {t is the amount of translation) in 
the Fourier domain. Hence translation ambiguity whenever present was corrected 
at the coarser levels only. 


3.7 Fourier coefficient truncation 

The Fourier transforms of general images have been observed to have most of their 
energy concentrated in a small region in the frequency domain near the origin and 
along both the dimensions of the frequencies. The reason for this can be that 
typically images have maximum region where the intensities change slowly. On the 
other hand sharp discontinuities like edges contribute to high as well as low frequency 
components. Since most of the signal energy is concentrated in a small frequency 
region, it should be possible for an image to be reconstructed without significant 
loss of quality from a fraction of the transform coefficients. Therefore to illustrate 
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the robustness of the reconstruction by any phase retrieval method, the Fourier 
coefficient truncation was done to 50% and to 12% of the Fourier modulus vzilues 
of the retrieved values. Rest all the transform coefficients were set to zero. It was 
hence shown that the reconstruction can be even attempted with fewer transform 
coefficients with a marginal trade off with quality and intelligibility. 
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Chapter 4 


Reconstruction experiments and 
results 


4.1 Fienup Reconstructions 


Test images taken up for phase retrieval by Fienup algorithm in conjunction with 
Hybrid Input-Output algorithm were two 128 x 128 pixel images 

• Planet Saturn 

• Figure 10 

The detailed program while is appended with this text, some salient aspects of the 
algorithm are 

• The algorithm was iterated for 11 cycles with each cycle comprising of 50 it- 
erations of Hybrid Input-Output algorithm and subsequent 10 iterations of 
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ORIGINAL 


INPUT GUESS 


CYCLE 1 





CYCLE 1 1 CYCLE 7 CYCLE 4 

Figure 4.1: Phase retrieval by Fienup algorithm. Clockwise two images (1) origi- 
nal Figure 10. (2) initial input guess, (3) to (6) are reconstructions at increasing 
iterations. 


4.2 Multiresolution wavelet based reconstructions 


Test images taken up for multiresolution decomposition and phase retrieval were 

Planet Saturn 
Figure 10 
Claire 
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(a) ORIGINAL 


(b) INITIAL GUESS 


(c) RECONSTRUCTION 


Figure 4.2. Phase retrieval by Fienup algorithm. (a) original Saturn (b) initial guess 
and (c) Final reconstructed image after 11 cycles 


Also since the decomposition to coarser resolution levels is wavelet based, following 
wavelets were chosen 

• Haar or Daubauchies-2 length( dbl) 

• Daubauchies 10 and 18 length (db5 and db9). 

• B-Spline Biorthogonal wavelets; BiorNd.Nr as Bior4.4(full length 10) and 
Bior6.8(full length 18). 

4.2.1 Characteristics of the wavelets used 

Orthogonal and compactly supported wavelets Daubauchies, Symlets and Coiflets 
form part of these.General properties of these class of wavelets are: 

• <f> and Ip exists and the analysis is orthogonal. 

• (p and Ip are compactly supported. 


NRMS Error 



Figtire 4.3: Reconstruction error for Saturn using Fienup algorithm for 11 cycles 

• ^ has a given number of vanishing moments. 

• Main nice properties of these wavelets are compact support and vanishing 
moments while difficulty is of poor regularity. 

• Daubauchies wavelets are characterised by their asymmetry. 

Biorthogonal and compactly supported wavelets B-Spline Biorthogonal wavelets 
form part of this class.General properties of this class are; 

• <f> function exists and the analysis is Biorthogonal. 
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NRMS Error 



Figure 4.4: Reconstruction error for Figure 10 using Fienup algorithm for 11 cycles 

• ^ and ijj both for decomposition and reconstruction are compactly sup- 
ported. 

• (f) and Ip for decomposition have vanishing moments. 

• 4> and Ip for reconstruction have known regularity. 

• Symmetry with FIR filters. 

• Desirable properties for decomposition and reconstruction are splitted. 

• Main difficulty is that the orthogonality is lost. 
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Figure Filrit ru»‘i}i( t»‘nt.s and fmjuency msponsps of selected Daubauchies 
wavelets. Cloekwi.M* frum top left (a) dbl ; High(*) and low pass(o) decomposi- 
tion wavelets (b) dh) |r) F:- ip!*-i;. y response of the wavelets; o shows dbl, * shows 
db5 and f sIhjws dh9 low pass tlc*€Oinposition wavelet frequency response, (d) db9. 

4.2.2 Algorithm 

Some saliimt aspot t.s uf the algorithm are 
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Figure 4.6: Filter coefficients and frequency responses of selected B-Spline Biorthog- 
onal wavelets. Clockwise from top left (a) Biorl.l ; High(*) and low pass(o) de- 
composition wavelets, (b) bior4.4 (c) Frequency response of the wavelets, o shows 
biorl.l, ★ shows bior4.4 and + shows bior6.8 low pass decomposition wavelet fre- 
quency response. (d) bior6.8. 

• Three levels of decomposition were carried out i.e.,m=0, . . . , 3 on all the three 
classes of images. 
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Figure 4.8; Compaxative Fourier moduli of bior4.4 (indicated on plot by o ) and 
bior6.8(*). All decomposition low pass characteristics. 

• The Fourier moduli at all the three levels of resolution were computed from 
the initial available Fourier modulus at zero(full) resolution. 
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Wavelet 

Mean error 

wavelet 

Fienup 

error 

Reduction 

in mean 

error 

Reduction 

in 

computations 


0.0295 

0.0324 

9% 

13% 


0.0277 

0.0324 

15% 

55% 

Db9 

0.0280 

0.0324 

14% 

>60% 

Bior44 

0.0296 

0.0324 

9% 

13% 

Bior68 

0.0266 

0.0324 

18% 

>60% 


Table 4.1: Reconstruction error comparisons <wavelet Vs Fienup> on image Saturn 

• The algorithm was started by taking a 16 x 16 array of random numbers for the 
images Saturn and Figure 10. For the Claire image the iterations were started 
by taking an input guess of 64 x 64 array of random numbers.The iterations 
were kept at 30 in all for the coarser resolutions (excluding full resolution for 
which 60 iterations for 11 cycles were taken). 

• Mask sizes were varied in a similar fashion as for the Fienup method i.e.,the 
initial mask was computed from the autocorrelation function at the full reso- 
lution. The support at subsequent resolutions are only the demagnifications 
appropriate to the resolution level. Although the decomposed levels should 
have support somewhat greater than half the support of the previous finer 
subspace, yet as an initial estimate, a tighter support constraint was taken 
which was subsequently loosened as the iterations progressed. 
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ORIGINAL 


INPUT 


( 1 ) 



( 4 ) ( 3 ) ( 2 ) 



( 5 ) ( 6 ) (y) 


Figure 1.9: Ph.uic retrieval by Wavelet algorithm using dbS wavelet.Clockwise from 
left tup (a) original Saturn (b) initial guess and (1) Final reconstructed image at res- 
olution 32x32 (2) reconstruction at level 64 x 64 resolution (3) to (7) reconstructions 
at increasing level of iterations for full resolution 

4.3 Ambiguity correction results 

As has been bought out in previous chapter, the reversal and translation are two 
ambiguities which are inherent in the phase reconstructed images when the only 
known input to start with is the Fourier modulus of the original image. The same 
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(a) ORIGINAL 


(b) 




(d) (c) 

Figure -I.IQ: Pluuse retrieval by Wavelet algorithm using dbl wavelet. Clockwise 
from left top (a) original Saturn (b) Final reconstructed image at resolution 32 x 32 
(c) reconstruction at level 64 x 64 resolution (d) reconstruction after 11 cycles at 
full resolution 

had been noticed in almost all the reconstructions whether partially or in full. In 
adrlition the translation ambiguity was found to be more prevalent amongst the 
two. These ambiguities were corrected in the manner as described in chapter 3. The 
results are illustrated in fig. 4.16. 
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(d) (c) 


Figun' -l-Il: Phsise retrieval by Wavelet algorithm using B-Spline Biorthogonal 
wavelet Bior4.4. Clockwise from left top (a) original Figure 10 (b) Final recon- 
structed image at resolution 32 x 32 (c) reconstruction at level 64 x 64 resolution 
(d) reconstruction after 11 cycles at full resolution 

4.4 Fourier coefficient truncation results 

Since most of the signal energy is concentrated in a small frequency region, efficiency 
of iterative transform and wavelet based reconstructions can be qualitatively mea- 
sured by truncation of their Fourier coefficients(i.e., limiting the images to certain 
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Figure 4.13: Reconstruction errors for Saturn using wavelet based approach for first 
seven cycles of iterations. (a) Using Bior6.8 wavelet (b) dbl (c) dbS (d) db9 and 
finally (e) Bior4.4. Convergence pattern remains same uptil 11 cycles. 


frequency regions only). To this effect the test images were taken and their Fourier 
coefficients were retained for upto 50% and 12%. The reconstruction error was al- 
though increased yet the image quality was preserved to a large extent of visual 
satisfaction. Quality and intelligibility were preserved without significant loss with 
a small fraction of the transform coefficients. Results of such comparisons are listed 
under table 4.1 for the Saturn image reconstruction with dbl wavelet and see figure 
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Figure 4.14: Reconstruction errors for Figure 10 using wavelet based approach for 
first .seven cycles of iterations, (a) Using Bior6.8 wavelet (b) dbl (c) dbS (d) db9 and 
finally (e) Bior4.4. Convsngf'nce pattern remains uncertain uptil 7 cycles wherein the 
nuisk was increased. Note that the input to the full resolution algorithm by dbl 
wavelet is lesser than all others. 

4.17 for the reconstructions as regards Saturn and figure 10. Error results in case 
of Figure 10 are similar in nature as of Saturn. Quality of reconstruction for all 
the wavelets in respect to Fourier reconstruction is satisfactory and the difference is 
not discernible by naked eyes due to artifacts of the viewing programs available and 
hence a representative figure only is being shown (figure 4.17). 



NRMS Error- 



Figure 4.15: R«construction errors for Figure 10 using wavelet based approach for 
8 ~ 11 cycles of iterations. (a) Using Bior6.8 wavelet (b) dbl (c) db5 (d) db9 and 
finally (e) Bior4.4. 

4.5 Discussions 

1. The iterative transform method of phase reconstruction is a robust method 
and will generally retrieve the phase from a wide variety of 2-D Fourier mod- 
uli. While most images can be reconstructed satisfactorily requiring only a 
few cycles of iterations, other more difficult objects can require many cycles 
without giving satisfactory convergence and reconstruction. Confidence that 
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Wavelet 

Mean error 

wavelet 

Fienup 

error 

Reduction 

in 

computations 

Dbl 

0.5600 

0.5590 

nil 

Dbo ! 

1 

0.5569 

0.5590 

67% 

Db9 ^ 

0.5548 

0.5590 

67% 


T;ibl{* 1.2: Reconstruction error comparisons < wavelet Vs Fienup> on Figure 10 



Figure 1.16: Ambiguity Corrections (a) Phase retrieved image obtained by using db5 
waveji't. Note the translation and reversal ambiguities, (b) Translation ambiguity 
eiurectetl image {<:) Reversal ambiguity correction superimposed to get the final 
pha,se retrii'ved image. 


the solution is the one and only true solution, can be increased by performing 
two or more trials of the algorithm, each time using different random numbers 
for the initial input. 

2. The Fienup algorithm has a limitation of not being able to fill up the details in 
an image during reconstruction. This seems to be due to the lack of adequate 
reconstruction in the high frequency regions resulting in blurring of edges when 
there is complete lack of phase data. Consequently the Fienup algorithm does 
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Tilhy 4.3; Reconstruction with truncated Fourier Coefficients 


Wavelet 

Percentage of 

retained 

coeffs 

N'ow error 

(previous error) 

Error 

Incr/Dcr 

Reconstruction 

I 


0.0407^0.0205) 

27.5% 

Satisfactory 

-({.> 

12% 

0.0831(0.0295) 

64.5% 

-do- 

Dh5 

30% 

0.0149(0.0277) 

-85.9% 

-do- 

-do- 

12*% 

0.0508(0.0277) 

51.2% 

-do- 

Dh9 

50% 

0.0498(0.0280) 

43.7% 

-do- 

-do 

12% 

0.0890(0.0280) 

68.5% 

-do- 

Biorli 

50% 

0.0471(0.0296) 

37.1% 

-do- 

-tlo- 

12% 

0.0852(0.0296). 

65.2% 

-do- 

BiorOH 

50% 

0.0411(0.0266) 

35.3% 

-do- 

i -do- 

j 

12% 

0.0863(0.0266) 

69.17% 

-do- 


not seetii to be suitable for detailed images where intelligibility measure is 
very exact ami mort; biased towards human visual perception. Human facial 
features are a <;<ise in point here. 

3. Tht‘ wavelet decomposition algorithm decomposes the search space into or- 
thogonal subspaces and hence is efficient for the following reasons: 

• It gives an initial input closer to the true solution as compared to the 
Fienup algorithm and hence chances of getting less residual error and 
adequate reconstruction are greater. 

• Avoids possible stagnation basins by reducing the low frequency error 
responsible for the slow convergence of the algorithm. 
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• Li'sser computational complexity while working in the subspaces and 
overall reduction in computations as compared to the Fienup algorithm. 

I. the wavelet decomposition algorithm presents an advantage over the the 

i'ieiiup algorithm, the choice of a particular wavelet depends on the following 

({iiaiity criterion; 

(a) Orthogonality of the bases. 

(b) Compactness of the wavelets. 

(<’} R.'gularily of scaling functions. 

(d) Number of vanishing moments. 

(e) Characteristics of the scaling function. 

1 he above factors influence reconstruction error , convergence and quality in 

the following manner: 

(a) Daubaehies class of wavelets are more compactly supported than the 
Biorthogonal spline wavelets. Hence the decompositions in the sub- 
spaces are more representative estimates than the Biorthogonal spline 
wavelets and therefore residual error will be less at lower resolutions for 
tlui Daubaehies class wavelets. 

(b) Wavelet support ( more) gives better and tighter support constraints as is 
obvious and hence in the same class of wavelets, a longer length wavelet 
will give lesser residual error. This is being borne by the results obtained 
wherein the reconstruction errors for lengthier filters (elements making up 
the support of the impulse response) are less. However although it is more 
likely that the initial error in case of longer length filter will be lesser as 
compared to shorter filter, yet the nature of iterative transform algorithm 
is such that the convergence pattern of the algorithm can not be clearly 
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ted tis the iterations progress. Many reconstruction experiments 
have l,een carried out and it can be concluded that although it is better 
to start with an input closer to the true solution, yet the nature of the 
' I can not be predicted by initial input alone. 

If) Regularity of Dauhachies(QA) class is poor when compared to Biorthog- 
onal spline wavelets(;V) and hence it is expected that it should introduce 
more discontinuities (artificial edges) at the coarser levels during decom- 
{jositioii. However since no PR is being attempted (no synthesis wavelet 
being n.sed) it is understandable that the lack of higher frequencies in the 
r«‘ct)nstructi«n is setting aside this advantage, the Biorthogonal spline 
wavelets would have had if the analysis-synthesis pairs were used for re- 
construction. 

(d) Linear phtise filters as compared to Non linear phase filters have better 
diH'omposition properties in so far as degradation of edges is considered. 
However the advantage Biorthogonal spline filters would have had in this 
respect has been again lost for the same reasons that only modulus in- 
formation of the low pass filtered image is being considered for phase 
retrieval. 

5. Ainhignity removal from the phase reconstruction algorithm has been at- 
tempted and it has been found that once the ambiguities are removed, they 
will not manifest again during the reconstruction iterations. This does not 
seem to be incidental as has been found by many experiments. The possible 
cause .seems to be the partially reconstructed phase characteristics. It has been 
experimentally established before [4] that ambiguities manifest themselves less 
for non symmetrically supported objects and are prominent in centrosymmet- 
ric objects because of strong symmetry around the center. Since during the 
course of iterations the objects are only partially reconstructed, it is likely that 
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!h«' nnihiginti«''s onc« corrected will not remanifest themselves. 
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(C) 


Figure 4.17: Reconstruction of images with truncated Fourier coefficients, (a) Phase 
retrieveni images with dbl wavelet (b) Images obtained by preserving 50% of Fourier 
cocfFici(>uts. All others are set to 0. (c) Same with 12% of the coefficients preserved. 
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Chapter 5 


Conclusions And Scope For Future 
Work 


Image reconstruction from the Fourier modulus data is an important problem and as 
yet not fully been solved satisfactorily. Iterative transform methods, since efficient, 
are the most sought after algorithms in this connection. The algorithms however 
suffer from major disadvantages : 

1. Lack of initial input guess. 

2. Slow convergence and stagnation of the algorithm. 

3. Intensive computational complexity. 

4. Reconstruction ambiguities. 

Any improvement to alter the influencing factors for image reconstruction 
error minimisation hence is being sought. In this the wavelet decomposition ap- 
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proach is being validated with the aim of suggesting an improvement in the initial 
input guess, the convergence rate and ambiguity removal. 

The approaches are being tested with various wavelet bases and on different 
classes of images to reduce inconsistency in measurement errors. The effect of the 
wavelets used for analysis was studied and results summarised. 


5.1 Scope for future work 


It is observed that the nature of wavelets used has an effect on resultant reconstruc- 
tions obtained by iterative transform algorithms. However it would be beneficial if 
the same is tested on a wider variety of wavelets to suggest an optimum wavelet for 
image analysis-synthesis for Phase retrieval problems. 

While estimation of Fourier modulus data was taken, it is assumed that 
the modulus is recorded in a noise free, blur free environment. The algorithm can 
be suitably modified to investigate the phase retrieval froip noisy Fourier modulus 
data. 


Since an improved algorithm (ambiguities have been taken care of) is at 
hand, it will be of interest to use both analysis and synthesis wavelet pairs (high and 
low pass) in a PR fashion to integrate the high frequency details in the reconstruction 
at the subspace levels. This approach may improve the image reconstruction error 
and more so the reconstruction quality. 
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